program dispersion_spectrum
  ! 理论频散谱计算程序
  ! 基于广义透射/反射系数方法
  ! 根据 libs/tds/src/MainProgram.f 改写
  !
  use constants
  use parameters
  use input_output
  use basic
  use grt_coefficients
  use matrix_operations
  use wavefield_calculation
  implicit none

  ! 局部变量声明
  integer :: i, j, ii
  real(kind=8), allocatable :: Vph(:)        ! 相速度数组
  complex(kind=8), allocatable :: omega(:)   ! 角频率数组
  complex(kind=8), allocatable :: zzz(:,:)   ! 结果矩阵
  complex(kind=8) :: integ(10)               ! 积分结果
  real(kind=8) :: kn 
  complex(kind=8) :: o                  ! 波数和角频率
  complex(kind=8) :: cs, cp                 ! 复速度
  integer :: file_unit = 11                 ! 输出文件单元号

  ! (0). 初始化基本常数（已由 constants 模块提供）
  ! (0). 初始化参数结构（已由 parameters 模块提供）
  call init_all_parameters()

  
  ! (1). 读取和检查输入参数
  call green_input()
  ! (2). 计算基本参数
  call green_basic()
  ! 输出基本参数
  print*, '!!!!!!!!!!!!!!!!!!! Basic Parameters !!!!!!!!!!!!!!!!!!'
  write(*,100) '              Time window width : Twin =', comp%Twin
  write(*,200) '     Sample points(time domain) : mt   =', comp%mt
  write(*,100) '       Source-receiver distance : R    =', source%r0
  write(*,100) '         Increment of time(sec) : dt   =', comp%dt
  write(*,100) '           Cutoff frequency(Hz) : fmax =', comp%fmax
  write(*,100) '     Increment of frequency(Hz) : df   =', comp%df
  write(*,100) 'The period structure length(km) : L    =', comp%L
  write(*,100) ' Increment of wavenumber(km^-1) : dk   =', comp%dk
  print*, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
  
100 format(1x,a,f12.4)
200 format(1x,a,i5)


  ! 分配数组内存
  allocate(omega(source%no))
  allocate(Vph(source%nv))
  allocate(zzz(source%no, source%nv))

  ! 构建频率和相速度数组
  do i = 1, source%no
    omega(i) = (source%mino + (i-1)*(source%maxo-source%mino)/(source%no-1.0d0)) * pi2
  end do

  do i = 1, source%nv
    Vph(i) = source%minv + (i-1)*(source%maxv-source%minv)/(source%nv-1.0d0)
  end do
  
  print *,lo, ls
  ! 主计算循环
  do i = 1, source%no
    
    ! 计算复速度
    do ii = 1, nly
      cs = (1.0d0 + log(omega(i)/(pi2*source%maxo))/(pi*media%qs(ii)) + aj/(2.0d0*media%qs(ii)))
      cp = (1.0d0 + log(omega(i)/(pi2*source%maxo))/(pi*media%qp(ii)) + aj/(2.0d0*media%qp(ii)))
      media%vs(ii) = media%vs0(ii) * cs
      media%vp(ii) = media%vp0(ii) * cp
    end do
    
    ! 对每个相速度计算响应
    do j = 1, source%nv
      kn = omega(i) / Vph(j)
      o = omega(i) - aj * comp%oi

      ! 计算位移响应
      call funval(kn, o, integ)
      ! ! 保存结果（使用瑞利波垂向位移分量）
      zzz(i, j) = integ(10)
      print *, i, j, zzz(i, j)
    end do
    stop

  end do


  ! 输出结果到文件
  ! 每行输出对应一个频率的所有相速度值
  open(file_unit, file=trim(Outname)//'.ds')
  do i = 1, source%no
    ! 将每个频率对应的所有相速度值写入一行
    write(file_unit, '(100e20.10)') (abs(zzz(i, j)), j = 1, source%nv)
  end do
  close(file_unit)

end program dispersion_spectrum